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Two phase galaxy formation: The evolutionary properties of galaxies 



o 
o 

> 

o 
m 

< 

o 

6 

(N 
> 
O 

m 

d 

T— I 

o 



M. Cooki'2*, E. Barausse^'\ C. Evoli\ A. Lapi^'^ ^ G.L. Granato*'^ 

^Astrophysics Sector, SISSA/ISAS, Via Beirut 2-4, 1-34014 Trieste, Italy 

^INAF, Osservatorio Astronomico di Padova, Vicolo dell' Osservatorio 5, 1-35122 Padova, Italy 
^Dept. of Physics, Univ. di Roma 'Tor Vergata', Via della Ricerca Scientifica 1, 1-00133 Rome, Italy 
^INAF, Osservatorio Astronomico di Trieste, Via G.B. Tiepolo 11, 1-34131 Trieste, Italy 
^Centre for Fundamental Physics, University of Maryland, College Park, MD 20742-41 1 1, USA 



23 November 2009 



ABSTRACT 

We use our model for the formation and evolution of galaxies within a two-phase galaxy 
formation scenario, showing that the high-redshift domain typically supports the growth of 
spheroidal systems, whereas at low redshifts the predominant baryonic growth mechanism is 
quiescent and may therefore support the growth of a disc structure. Under this framework we 
investigate the evolving galaxy population by comparing key observations at both low and 
high-redshifts, finding generally good agreement. By analysing the evolutionary properties of 
this model, we are able to recreate several features of the evolving galaxy population with 
redshift, naturally reproducing number counts of massive star-forming galaxies at high red- 
shifts, along with the galaxy scaling relations, star formation rate density and evolution of the 
stellar mass function. Building upon these encouraging agreements, we make model predic- 
tions that can be tested by future observations. In particular, we present the expected evolution 
to z=2 of the super-massive black hole mass function, and we show that the gas fraction in 
galaxies should decrease with increasing redshift in a mass, with more and more evolution 
going to higher and higher masses. Also, the characteristic transition mass from disc to bulge 
dominated system should decrease with increasing redshift. 
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1 INTRODUCTION 

The current paradigm of cosmological structure evolution is out- 
lined by the ACDM model: Providing a remarkably successful 
framework for interpreting a wealth of observations of cosmic 
structure evolution over the majority of the duration of the Uni- 
verse. This model is capable of reproducing the cosmic microwave 
background radiation fluctuations (Spergel et al. 2007), the large 
scale clustering of galaxies (Eisenstein et al. 2005, & references 
therein), the cosmic shear field measured through weak gravita- 
tional lensing (Hoekstra et al. 2006 & references therein), small 
scale power spectrum of Lyman-alpha forest sources (Jena et al. 
2005), the properties of galaxy clusters (Allen et al. 2004 & ref- 
erences therein) along with several other key observations of large 
scale cosmological structures. However, despite these merits, on 
galaxy scales the assembly of baryonic material within virialised 
dark matter (DM) haloes has had more mixed successes: Due to 
the complex processes, often non-linear and dissipative, which op- 
erate on scales well below the resolution of the model ('sub-grid' 
physics), in order to model the evolution of baryonic material 
within DM haloes one is required to adopt analytic prescriptions 
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and make several important assumptions concerning the geometry 
of the forming system (see Zavala, Okamoto & Frenk, 2008). 

Early endeavors to model the cosmological evolution of lu- 
minous structures came from White & Rees (1978) and Blumen- 
thal et al. (1984), whereby galaxies form when gas cools and con- 
denses within the centres of hierarchically evolving DM haloes. At- 
tempts to model and interpret the evolutionary properties of galax- 
ies within the first generations of semi analytical models (SAMs) 
showed promising qualitative agreements to observations (Kauff- 
mann et al. 1993, 1998, Cole et al. 1994, 2000, Somerville & Pri- 
mack, 1999). However, in the past decade it has become clear be- 
yond reasonable doubt that significant tensions between SAM pre- 
dictions and fundamental observations exist, most notably in three 
major areas: Firstly the issue of 'overcooling' ('quenching'), which 
has several manifestations, large DM haloes are observed to be low 
in baryonic mass and contain typicaly 'red and dead' early-type 
galaxies, resulting in a sharp cutoff in the high-mass end of the 
stellar mass function, unlike that for the DM haloes (Bell et al. 
2003a, Benson et al. 2003, see Somerville et al. 2008b for a dis- 
cussion). Secondly, 'downsizing', or 'anti-hierarchical' evolution 
of baryonic structures (Cowie et al. 1996), whereby massive star 
forming systems and associated SMBHs shined mostly at high red- 
shifts, while smaller objects show longer lasting activity (see also 
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Figure 1. The evolution of the galaxy halo mass function (GHMF) shown 
with the solid curves, from < z < 5 and the Sheth-Tormen mass function 
(STMF) shown with the dotted curves. As can be seen, in all redshift ranges 
the GHMF resembles the STMF until approximately 10^^ where the 
probability for multiple galaxy halo occupation grows and thus the single 
galaxy mass probability diminishes, resulting in the exponential cutoff. 

Fontanot et al. 2009 for details) which appears contrary to naive ex- 
pectations for the 'bottom up' growth of DM structure. Finally, the 
'dwarf galaxy', or 'substructure' problem, whereby the number of 
low mass galaxies predicted by models is significantly more than is 
observed (see Mo et al. 2005, Moore et al. 1999). 

Theoretical attempts to interpret these somewhat puzzling 
properties of galaxies motivated a second generation of SAMs, 
which evoked strong feedback from a central supermassive black 
hole (SMBH) in order to quench star formation at late times by 
suppression of cooling, predominantly in the larger galaxies (see 
Croton et al. 2006, Bower et al. 2006, Baugh et al. 2006, see also 
Granato et al. 2004), generating a marked improvement over pre- 
vious incarnations, but several tensions remained (see Monaco et 
al. 2007). The current state-of-the-art SAMs include the energetic 
effects of growing central SMBH, the effects of hot and cold accre- 
tion (Dekel & Birnboim, 2006, Cattaneo et al. 2006, Somerville 
et al. 2008b) or a flat stellar initial mass function (IMF) during 
starburst activity (Baugh et al. 2005), the suppression of cooling 
and collapse due to an ionizing UV background (see Gnedin, 2000, 
Somerville 2002, Benson et al. 2002), thus steadily increasing the 
degrees of freedom in order to improve agreement with observa- 
tional constraints. Progress is currently being made in developing 
SAMs with added layers of physical descriptions, including spa- 
tially resolved modeling (see Stringer & Benson, 2007, Dutton & 
van den Bosch, 2009, Cook et al. 2009a (hereafter C09a)), multi- 
phase ISM physics (Dutton & van den Bosch 2009, Cook et al. 
2009b (hereafter C09b)) in order to increase predictability of mod- 
els without significantly increasing their number of free parameters. 

Since the first generation of SAMs were developed, observa- 
tional studies have undergone many revolutions due to increased 
sensitivity, increased wavelength coverage, and automated survey 
methods. Many of the observational constraints coming as a sur- 
prise to the community: At low redshifts, detailed constraints on the 
stellar mass function (Cole et al. 2001, Bell et al. 2003a) improved 
model parameter refinements, however, analysis of high redshift 
star-forming galaxies {z ^ 2) opened a window to study the prop- 
erties of galaxies when the Universe was under 20% it's current 
age, lyman break galaxies (Steidel et al. 1996), Sub-mm galaxies 



(Small et al. 1997) and Ly-a galaxies (Hu, Cowie & McMahon, 
1998) were generally interpreted as being dusty starbursting sys- 
tems, with detailed analysis showing that the star formation rate 
density of the Universe at 2: > 2 remains flat (in contrast with the 
original determination of Madau et al. 1996). Also, measurements 
of the mass distribution of high-z galaxies revealed a substantial 
population of extremely massive galaxies at z > 1 (Cimatti et al. 
2002, Drory et al. 2003, Kodama et al. 2004, Bundy et al. 2005) 
in sharp contrast to the original hierarchical picture of structure 
growth. Current observations of stellar mass now extend to 2 ~ 5 
(Drory et al. 2005, Fontana et al. 2006, Eisner et al. 2008, Perz- 
Gonzlez et al. 2008, Marchesini et al. 2009), and theoretical models 
must attempt to interpret these results physically whilst simultane- 
ously making predictions about the black hole growth (Hopkins et 
al. 2006) for which observations are complete to high redshifts, and 
the scaling relations of galaxies and SMBHs (see Woo et al. 2008) 
along with the 'archeological' constraints on the evolutionary prop- 
erties of galaxies (see Gallazzi et al. 2005). It has been shown that 
theoretical models have had mixed successes, with no model cur- 
rently able to consistently predict all observations (see Kitzbichler 
& White, 2007, Marchesini & van Dokkum, 2007, Somerville et al. 
2008b, Fontanot et al. 2009). 

Despite several differences in the detailed 'sub-grid' recipes 
adopted by different groups, current SAMs all follow the same gen- 
eral framework as originally proposed by White & Rees, 1978 and 
adopt the same original assumptions: (i) Gas cooling and conden- 
sation with in DM haloes, at any epoch, results in the dissipation- 
less formation of a self-gravitating gaseous disc which undergoes 
mild star formation, (ii) The main driver for starburst activity is the 
merging of these gas-rich discs (wet mergers) which also provides 
the main channel for the formation of spheroidal structures (Cole 
et al. 1991). The resultant 'disc merger' framework provides the 
basis for most current SAMs (see Somerville et al. 2008b for a re- 
view). However, in our view these strict assumptions may be the 
underlying cause of several tensions between models and observa- 
tions (see Mo & Mao, 2004, C09a), notably the tendency for bary- 
onic material to follow the hierarchical evolution of DM haloes, the 
difficulty in producing massive galaxies at early times which later 
passively evolve, and archeological issues relating to the structure 
of DM haloes and the observed baryon fraction in galaxies (see the 
aforementioned references). 

Within this work, we develop the model outlined in C09a, 
C09b, where, motivated by the above-mentioned tensions between 
theoretical models and observations we proposed a model which 
differs substantially from the standard 'disc merger' framework: 
We envisage that the fundamental dichotomy between galactic 
spheroid and disc components is a manifestation of two distinct 
modes of the evolution of baryonic matter, ultimately driven by 
the two-phase structural evolution of DM haloes (see Zhao et al. 
2003a, Mo & Mao, 2004, Diemand et al. 2007, C09a, C09b): An 
early 'fast collapse'' phase, where the DM core structure is con- 
structed through a series of violent merger events, corresponding 
to an epoch where baryonic material effectively dissipates angu- 
lar momentum upon collapse to directly form a spheroid-SMBH 
system, and a late 'slow collapse'' phase, where potentially large 
amounts of material are added to the halo outskirts little affecting 
the central regions, giving rise to the quiescent growth of disc struc- 
tures around the pre-formed spheroids. 

Dark matter mergertrees outline the merging rates of DM 
haloes and are well constrained by simulations. However, ultra 
high-resolution simulations are required in order to analyse the 
structure and substructure evolution within the merging DM haloes. 
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Figure 2. The stellar mass function evolution: Locally (top left) plotted against the determinations of Cole et al. 2001 showing good agreement with both the 
cutoff and normalisation. At z > compared to Marchesini et al. 2009, showing the z » 1.5 relation (top right), reproducing the high mass cutoff and knee, 
but slightly overproducing the number of low-mass galaxies, the z Ri 2.5 relation (bottom left), showing an overall good fit within the observational range, and 
the 2 3.5 relation (bottom right), showing good agreement within the data ranges, but not showing a good approximation to a Schechter function fit. Yellow 
shaded regions represent the errors in the functional fits (orange lines), and the blue error bars in the model outputs representing the Poisson uncertainty in the 
mean averages. We note that slight discrepancies in the low mass end may be somewhat attributed to observational biases, and environmental effects which 
we neglect to model (see §4 for a further discussion.) 



Thus, until recently, oversimplified analytical recipes are com- 
monly used (Chandraseakhar, 1943), not accounting for several im- 
portant effects. Recently, increased numerical resolution and sub- 
structure analysis has allowed for some advances in determining 
the evolution of subhaloes after they have entered a parent halo 
(which is of upmost importance for baryonic physics), showing that 
in general, the evolution of the structure of a galaxy-sized DM halo 
evolves in two-phases. 

More specifically, analysis of the cosmological evolution of 
virialised structures is long standing, from observational clustering 
studies, through detailed cosmological simulations (Springel et al. 
2005), and monte-carlo algorithms tuned to reproduce these results 
(see Parkinson et al. 2008 & references therein), however, until 
relatively recently, determining the detailed evolution of substruc- 
ture within DM haloes after they merge has been somewhat over- 
looked. Recent increases in numerical resolution within N-body 
simulations have begun to analyse the detailed structural evolution 
of haloes within cosmological volumes (Zhao et al. 2003a, 2003b, 



Diemand, Kulhen & Madau, 2007, Hoffmann et al. 2007, Ascasibar 
& Gottloeber, 2008), showing that two distinct phases of structural 
evolution are found, an early 'fast collapse' phase, followed by a 
late 'slow collapse' phase. This has also prompted several works to 
show how typical double power-law DM halo density profiles may 
be generated (see Lu et al. 2006, Lapi & Cavaliere, 2009). This the- 
oretical idea has also been hinted upon in the Millennium Galaxy 
Catalogue bulge-disc decomposition analysis of Driver et al. 2006. 

Motivated by these issues, within this contribution we expand 
the model presented in C09a and C09b, which comprises a natu- 
ral extention to the spheroid-SMBH co-evolution model presented 
in Granato et al. 2004, (see also Granato et al. 2001, Lapi et al. 
2006, 2008, Mao et al. 2007), and focus on several 'problem plots' 
for current SAMs under the disc-merger framework. We essentially 
inherit the results of these papers here. By self-consistently out- 
putting galaxy properties at several different redshits for a repre- 
sentative sample of galaxies generated by our model, we are able 
to model the evolutionary development of baryonic material within 
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Figure 3. The number density of galaxies with stellar masses M > 
10^^ Mq and M > 10^^ Mq produced by Drory et al. 2005 and com- 
pared with model prescriptions. Yellow shaded regions represent the error 
in the functional fits (orange lines), and the blue en'or bars in the model out- 
puts representing the Poisson uncertainty in the mean averages. Due to the 
large star-burst activity at high-z we are able to construct massive galaxies 
readily and produce good agreements to the data. 



our models, comparing the fundamental relations in order to con- 
strain the key physical mechanisms governing galaxy formation. 
Focusing on the evolution of the galaxy scaling relations for discs 
and spheroids, the evolution of the mass functions for both SMBHs 
and galaxies, the cosmological star formation rate density, the cos- 
mological evolution of the most massive galaxies and the archeo- 
logical stellar populations of local galaxies. 

The plan of this paper is as follows; in §2 we overview the 
physical model, highlighting important points and modifications 
to previous works, in §3 we describe the methods in order to ex- 
tract observable quantities and present the results for the evolv- 
ing galaxy population, we conclude and summarise our findings 
in §4, highlighting the successes and limitations of our approach. 
Throughout the paper we adopt the standard ACDM concordance 
cosmology, as constrained by WMAP 5-year data (Spergel et al. 
2007). Specifically, we adopt a flat cosmology with density param- 
eters Q,M = 0.27 and JIa = 0.73, and a Hubble constant Ho = 70 
km s~^ Mpc~^. 



2 OVERVIEW OF THE MODEL 

We refer the reader to C09a and C09b for a detailed description of 
the model details. However, in order to preserve clarity we review 
the main model features here. 

For the dark matter merging and accretion evolution, we adopt 
an extended Press-Schechter formalism based on the binary merg- 
ertree of Cole et al. 2000, as modified by Parkinson et al. 2008. This 
algorithm has been shown to reproduce halo merging and accre- 
tion statistics obtained from cosmological numerical simulations 
(Springel et al. 2005). We use these mergertrees by extracting the 
main-progenitor mass accretion history (MAH) beginning at z = 
and moving to progressively higher redshifts in the mergertree, tak- 
ing the largest progenitor branch at each merger event. 

It has been shown in high-resolution simulations, and inves- 
tigated in analytical dark matter studies, that the structure of DM 
haloes evolves in two phases, a 'fast' accretion phase at high z 



Figure 4. The cosmological star formation rate density evolution. Model 
comparisons with the Somerville et al. 2001 compilation, showing good 
agreements to observations between < 2; < 5, also plotted are the con- 
tributions from the spheroid (red) and disc (blue) components. The yellow 
shaded region represents the errors in the functional fits (orange line), and 
the blue error bars in the model outputs representing the uncertainty in the 
mean averages. These results show, in broad agreement with observations, 
that elliptical galaxies (and galaxy bulges) form early, and discs form late, 
with a typical transition at z RJ 1. 



and a 'slow' accretion phase at low z. These two phase are re- 
flected in the redshift evolution of the concentration parameter 
c(z), which characterises the halo structure. In our work, we cal- 
culate c{z) by means of recent simulation results for the z = Q 
mass-concentration relation (Maccio et al. 2007), coupled to the 
evolutionary evolution reported in Zhao et al. 2003a, 

[ln(l + c) - c/(l + c)]c-'" cx H{zf''K'Ui,{zf-°' , (1) 

where a is a piecewise function (Zhao et al. 2003a). The knowl- 
edge of the evolution of c{z) then allows us to distinguish the slow 
and fast DM accretion phases, which we associate with two growth 
mechanisms for the baryonic sector: the fast DM accretion phase 
giving rise to the the formation of bulges, and the slow DM accre- 
tion phase giving rise to the the formation of discs (Mo & Mao, 
2004, C09). 

More specifically, in order to model the baryonic evolution, 
we start at a redshift at which the virial mass (given by the MAH) 
reaches the so-called cooling mass, i.e. the virial mass at which 
Tvir ~ W^K. In fact, below this temperature the Lya cooling 
becomes inefficient and baryonic structures cannot form. Follow- 
ing this, over each redshift increment, we allow hot gas to accrete 
onto the DM halo with rate Mini = /coiiAfvir, where /coii is the 
baryonic collapse fraction in the presence of an ionizing UV back- 
ground (Gnedin et al., 2004, Somerville et al., 2008b) 

m 



/coll(A'/vir, z) 



(2) 



(l + 0.26M/(z)/Afvir)=' 

(Mf{z) is the filtering mass at a given redshift: see Kravtsov, 
Gnedin & Kyplin, 2004, Appendix B). Also, we include the effects 
of cold accretion flows, shown to be the predominant mechanism 
leading to the formation of low-mass systems. Below a critical mass 

Mc = Ms max[l, lo^-^t^-^^)] , (3) 

where Ms — 2 x 10^^ Mq and Zc — 3.2, we assume that all gas 
accreted onto DM halos is not shock heated to the virial tempera- 
ture of the DM halo, but streams in on a dynamical time (see Dekel 
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Table 1. Values of the free parameters of our model. 



Description 


Symbol 


Fiducial value 


Impact on this work 


SN feedback efficiency (bulge) 


<:SN,b 


0.5 


Strong 


SN feedbacic efficiency (disc) 


CSN.d 


0.8 


Strong 


Reservoir growth rate 


^res 


10"'"'Ai?,;,yr-l 


Strong 


QSO feedbaclc efficiency 


/h.QSO 


10-4 


Strong 


Radio feedbacif efficiency 


/h, radio 


1 


Strong 


Viscous accretion rate 


kacc 


10-2 


Weak 


Radio mode accretion rate 


^radio 


6 X IQ-^Moyr-l 


Weak 



& Birnboim, 2006, Dekel et al., 2009, Cattaneo et al, 2006). Thus, 
in halos below this mass the collapse happens on the dynamical 
timescale of the system (tcoii = tdyn), whereas in halos above this 
mass tcou ~ riiax[idyn! ^cooi], where the cooling timescale tcooi is 
computed in a standard way, assuming material is shock heated to 
the virial temperature. The effects of this cold accretion is to en- 
hance star formation at high redshifts relative to the scenario where 
all material is shock heated. 

In order to model the baryonic evolution, we suppose that the 
hot gas phase collapse gives rise, as we have already mentioned, to 
bulges during the fast accretion phase and to discs during the slow 
accretion phase: 



Afcoll 



A^b,gas 
-^^d.gas 



lz> Zt] 
[z<zt] 



Zt being the transition redshift between the slow and fast accretion 
phases. Thus, we naturally output the growth of a spheroid structure 
followed by the growth of a disc structure around the pre-formed 
spheroids. 

For z > Zt, as gas collapses into the bulge, bursts of star for- 
mation occur which force, by radiation drag, part of the cold gas 
onto a circumnuclear reservoir with low angular momentum, at a 
rate 



(4) 



The cold gas in this reservoir then becomes eligible to feed a central 
seed supermassive black hole at an accretion rate 



Mbh.Qso = min[Mvisc, Medd] , 



(5) 



where A/, rid is the Eddington rate and the viscous accretion rate is 
parameterized as 



■ ^acc 



G V Afbh 



(6) 



where fcacc 10~^ is a free parameter with little impact on our 
results. Feedback on the growth of baryonic structure comes from 
two processes. On the one hand, supernova (SN) explosions trans- 
fer significant energy into the cold ISM, causing it to be re-heated 
and ejected from the system. Therefore, by considering energy bal- 
ance in the ISM, we assume that supernova feedback is able to re- 
move gas from the bulge with efficiency esN.b (ranging from to 
1, with esN.b = 1 meaning that all of the SN explosion energy is 
adsorbed by the ISM). This mechanism is most effective in the low 
mass systems, which presents shallow potential wells from which 
the ISM can easily escape due to SN explosions. 

On the other hand, the QSO activity of the central SMBH 
ejects hot gas and bulge cold gas from the system with an efficiency 



/ii.QSO (/ii,QSO ranging from to 1): 

2 Lh Afb.gas 



M 



QSO 

bjgas 



= /h.QSO 



3 Mhot + A^b.gas 



M, 



QSO 



Jh.QSO- 



Mv,, 



3 cr2 A^hot + A^b.f 



(7) 



(8) 



where a — 0.65Kir (Kir being the halo virial velocity), while 
Afb.gas and Afhot are the masses of the gaseous bulge and of the hot 
gas phase. This effect is most effective in the large mass systems, 
where QSO activity is strong. 

In addition to the QSO accretion channel, we assume, follow- 
ing Croton et al. 2006, that the SMBH accretes mass also through 
a quiescent "radio-mode" at a rate 



M, 



bh radio 



adi 



/_Mbh_ 



(9) 



\10^MqJ \0.1 J V200km/s 

where /hot is the halo mass in the form of hot gas and fc^adio = 
6 X IQ-^Afoyr"^ as in Croton et al. 2006. Because of the small 
value of fcradio, this mode does not contribute significantly to the 
SMBH mass evolution. However, following Croton et al. 2006, we 
assume that the efficiency /ii, radio with which the energy emitted 
by the SMBH in this mode is adsorbed by the hot gas phase is 
exactly 1 (i.e., all of the radio mode emission is adsorbed by the 
hot gas phase). 

As the DM halo enters into the relatively quiescent 'slow' ac- 
cretion phase at z < zt, the DM halo core potential becomes sta- 
bilised and we suppose that conditions become sufficient to support 
the growth of a disc through dissipationless collapse. Thus, gas en- 
tering into the DM halo conserves angular momentum and joins a 
gaseous disc, for which we assume an exponential surface density 
profile with scale radius is calculated following (Mo, Mao & White, 
1998, equation 29, and C09b equation 31). We adopt for simplicity 
the same for both the gas and stellar discs, but we have tested 
that a somewhat larger scale-length for the gas rd.gas = 'i--5rr,stars 
(e.g. Somerville et al 2008), does not yield any significant differ- 
ence in the results discussed here. Star formation in these gaseous 
discs is expected to take place in molecular clouds (see C09b, sec- 
tion 3.3 for more details on the star formation law that we use) and 
gives rise to a stellar disc, for which we assume an exponential sur- 
face density profile with the same scale radius as the gaseous disc. 

It is known that when discs become self-gravitating they are 
likely to develop bar instabilities, get disrupted and transfer ma- 
terial to the spheroidal component (Christodoulou, Shlosman & 
Tohline, 1995). We therefore assume that a stellar or gaseous disk 
is stable if 

K(2.2rd) 



: stars , gas , 



(10) 



where af.ff = 1.1 and a^^^ = 0.9 [see (Mo, Mao & White, 1998) 
and references therein]. If we find that discs become unstable, we 
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Figure 5. The evolution of the stellar mass Tully-Fisher relation, showing model comparisons with Bell & de Jong, 2001 for local galaxies (left panel) and 
finding excellent agreement across the observational range, with Conselice et al. 2005 (center panel), finding again good agreement and little evolution from 
the 2 = relation, and finally a tentative comparison with the 2 = 2 results of Cresci et al. 2009 (right panel), finding an offset, see §3.4 for details. Yellow 
shaded regions represent the scatter in the relation, containing the 95th percentiles, and the mean values are given by the orange lines. 



assume they get disrupted in a dynamical time and transfer their 
material (either stars or gas) to the bulge components. 

Feedback on the disc growth comes again in two fashions. On 
the one hand, analogously to the bulge case, we assume that super- 
nova explosions can remove gas from the disc with efficiency eg N,d 
(ranging from to 1). Again, this mechanism is only efficient for 
small systems. On the other hand, QSO activity is not generally 
present in the slow accretion phase, unless the gaseous disc frag- 
ments due to bar instability into a spheroidal gaseous component, 
which immediately forms stars, feeding the reservoir and, through 
it, the SMBH. However, the radio mode feedback is still present 
and removes hot gas from the system, thus quenching the collapse 
of the hot gas into the disc cold gas and indirectly suppressing disc 
star formatiorQ 

Finally, in order to account for adiabatic halo response, we 
take the standard prescription of Blumenthal (1986). In particular, 
denoting by Mx{r) the mass of a given component 'X' (X — h 
for the bulge, X = d for the disk, X = DM for the dark matter 
and X — res for the reservoir) enclosed by a radius r, from the 
angular momentum conservation one obtains 

M^{n)r^ = Mf{r;)r; , (11) 

where Vi and r/ are respectively the initial and final radius of the 
shell under consideration, the initial mass distribution Mi[ri) is 
simply given by the NFW density profile, while Mf{rf) is the final 
mass distribution. Also, mass conservation easily gives 

Mjirf) = Afd(r/) + A'/b(r/) + Mj^mirf) + M^csirf) = 
Md(r/) + Mb(r/) + M,es(r/) + (1 - fgai)M,{r,) , (12) 

where fgai = Mgai/M^i, (with Mgai = Md + Mb + M^cs)- By 
assuming spherical collapse without shell crossing, one can adopt 
the ansatz r/ = Fri, with F = const (Blumenthal, 1986), and 
Eqs. nil and l ll2t can be solved numerically for the contraction 

^ In our model we assume the QSO and radio mode feedback do not re- 
move cold gas from the disc, due to the small geometric cross section of the 
disc relative to the SMBH emission. 



factor r. However, in order to be able to mitigate or even switch 
off the halo adiabatic contraction, we modify by hand the relation 
between and r/ and assume, as in Dutton et al. 2007 

r/^F^n, (13) 

where /i is a free phenomenological parameter. Therefore, /i = 1 
corresponds to adiabatic contraction as in Blumenthal, 1986, while 
^ — completely switches off adiabatic contraction. 



3 RESULTS 

In order to make comparisons with observations, we produce a sta- 
tistical sample of approximately 1000 galaxies with z = virial 
masses in logarithmic increments in the range 9.5 < log{Mv) < 
13.5. In order to account for the cosmological abundances of galax- 
ies, we assign each DM halo a weight using the galaxy halo mass 
function (GHMF). This was originally derived by Shankar et al. 
2006 in order to account for the one-to-one relationship between 
galaxies and their host DM haloes. Essentially it is derived using 
numerically and constrained DM halo mass functions (see Sheth 
& Tormen, 2002 (STMF), Jenkins et al. 2001) but accounting for 
the halo occupation distribution (HOD) within DM haloes, which 
is unity for the majority of galaxy hosting systems, but rapidly in- 
creases in haloes within haloes with A/„ > 10^'^ Mq. Thus, to ac- 
count for this, we subtract the subhalo mass function (SHMF) as 
derived in van den Bosch et al. 2005. The benefit of using these 
parameterisations is that they are also defined at 2: > and thus 
may be used to extract the number densities of galaxies at higher 
redshifts, within this work, we use the following: 

™^(^'^)-^f(l^(|)'"-p(-|) '''' 

Where a, /3 and 7 are given in van den Bosch et al. 
2005, and — ra/M are the normalised subhalo masses. Thus 
GHMF(M,z) = STMF(M,z) - SHMF(M,z) is essentially 
identical to the STMF at masses below 10^"^ at 2; = 0, with an ex- 
ponential cutoff at higher masses due to the dominance of groups 
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Figure 6. The disc size as a function of stellar mass at z = (left panel) and z = 1 (right panel) as compai'ed to the results deiived in Somerville et al. 
2008a and based upon SDSS galaxies (Blanton et al., 2005) locally, and GEMS galaxies (Rix et al., 2004, Caldwell et al., 2005) at high redshifts. Finding 
good agreements between model and observation at z = but noting a slight offset when compared to the z = 1 sample. Yellow shaded regions represent the 
scatter in disc scale lengths at fixed stellar mass, containing the 95th percentiles. Orange lines show the mean values for the size-mass relation. 



and clusters of galaxies. In Fig. 1, we show our derived GHMF at 
different redshifts, as can be seen, at 2: = the GHMF is identical 
to the SFMF, but at A/„ > IO^Mq the GHMF exponentially 
drops, having a negligible probability at Mv > 10^^ Mq. Since 
clusters form at relatively late times within the standard hierarchi- 
cal picture, we do not see significant evolution in this cutoff mass to 
high redshifts, however we do see the typical evolution in the mass 
function. 

3.1 Galaxy stellar mass function Evolution 

One of the fundamental constraints on the physical mechanisms 
governing the evolution of luminous matter in galaxies is encoded 
within the stellar mass function, since its shape holds an imprint of 
the underlying physics which dominates on different mass scales. 
Typically, the mass function is fit accurately by a 'Schechter' func- 
tion (Schechter, 1976) with a low mass power law slope a, a char- 
acteristic mass M,t, and normalisation $*. 

It is generally understood that the low mass power law slope 
may be matched with a combination of ionizing UV background 
suppression of infall and supemovae feedback, since the poten- 
tial wells of their host DM haloes are relatively shallow and can- 
not capture and retain baryonic material (see Benson et al. 2002). 
Whereas the bright end has proved to be more of a challenge, and 
is now understood as a combination of cooling inefficiencies cou- 
pled with multiple occupation and strong energetic feedback from 
a central SMBH (Granato et al. 2004, Bower et al. 2006, Croton et 
al. 2006). These theoretical predictions, however, have been shown 
to show some discrepancies at higher redshifts (see Marchesini et 
al. 2008, Fontanot et al. 2009, Kitzbichler & White, 2007, De Lucia 
& Blaizot, 2007) 

From an observational perspective Cole et al. 2001, and Bell 
et al. 2003a used near-IR colours in order to determine the stel- 
lar masses, however, more recent approaches model masses using 
multi-wavelength approaches (Drory et al. 2004, 2005, Fontana et 
al. 2006, Perez-Gonzalez et al. 2008, Marchescini et al. 2008), ex- 
ploiting broad-band photometry to compare with libraries of syn- 
thetic spectral energy distributions (SEDs) which output the best 



fitting photometric redshift, stellar mass and SFR. Thus, the deter- 
minations of stellar masses is subject to several model-dependent 
uncertainties and simplifications (such as a smooth star formation 
history interspersed stochastically with starburst events, unlike the- 
oretical models, which typically exhibit complex histories) and the 
results are therefore subject to several potential biases (see March- 
esini et al. 2008 for an extensive analysis). 

Within this work, we compare model predictions between 
< z < 4 with the results of Cole et al. 2001 for local galaxies, 
Marchesini et al. 2008 for (z > 0) populations as shown in Fig.2. 
Locally we find a good agreement with the observed mass function 
in both high mass cutoff and normalisation, consistently reproduc- 
ing the observations down to Mstars ~ 10^" A/0 . However, in the 
lowest mass systems, we do find a slight discrepancy, overproduc- 
ing the number of low-mass galaxies (an effect which may also be 
seen in other SAMs, see Fontanot et al. 2009 Fig.l). At higher red- 
shifts we are able to generate a close match to the high mass cutoff 
up to 2: ~ 4, we view this as a notable success of our model since 
several other current SAMs find this difficult (see Marchesini et al. 
2009, Fig. 13), typically under-predicting the cutoff mass and over- 
producing the number of low mass galaxies. We do, however, find 
that at 2 > our model generates too many low mass galaxies 
which manifests clearly in the lowest mass systems, however, these 
mass scales are beyond the range of the observational constraints 
and thus it remains unclear as to the true faint-end slope at higher 
redshifts. 

Qualitatively, we may view the successful reproduction of 
the high-mass cutoff as a manifestation of the direct formation 
of spheroid-SMBH systems at high redshifts. Very early collapse 
onto spheroid structures without prior disc-formation results in the 
growth of large galaxies at early times, allowing for the high-mass 
end of the mass function to be in-place already at 2 > 4, in broad 
agreement also with the concept of 'cosmic downsizing', however, 
we do find that the overall Schecter function fit poorly describes the 
model and observation at z = 3.5, and therefore a comprehensive 
analysis of the exponential cutoff cannot yet be achieved. 
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Figure 7. The spheroid black-hole scaling relations locally and at z > 0. The local M-sigma relation as output by model and compared to the results of 
Tremaine et al. 2002, Greene et al. 2004 (top left) shows agreement with observations over the entire constrained range, and at z = 0.5 (top right) we 
compare model predictions with Treu et al. 2007, Woo et al. 2008, finding no significant evolution. For the mass relationship between bulge and SMBH, we 
compare at 2 = to the observations of Haring & Rix, 2004, finding promising agreements (bottom left), and we predict the form of this relation at z = 1, 
also comparing to the z = sample, and finding no significant evolution (bottom right). Yellow shaded regions represent the errors in the functional fits to 
observations (orange lines) in all panels. 



3.2 Massive galaxy number count evolution 

In order to further quantify the growth and evolution of the largest 
galaxies in the Universe, we compare model predictions with the 
massive galaxy number density evolution observations between 
< z < 5 by Drory et al. 2005, who found that the number den- 
sity of the most massive systems evolves in a manner similar to 
he evolution of lower mass systems and are present at all redshifts 
within their range. They derive this result by obtaining the stel- 
lar mass function for a sample of multicolour observations in the 
FORS Deep Field (Heidt et al. 2003) and the GOODS-south sur- 
vey (Giavalisco et al. 2004). By fitting Schechter-functions to the 
observations at a number of redshift intervals and integrating the 
results they determine the total stellar mass density evolution, and 
the galaxy number count evolution. Despite the significant uncer- 
tainties due to model-dependent stellar mass determinations, and 
Schecter-fitting, a striking relation has been obtained, showing that 
the largest mass systems are being formed at all redshifts within 
their range, and at « = 5 a significant number of large mass sys- 
tems are already formed. 



Outputting model predictions, we see in Fig.3 that we reliably 
reproduce the number densities of the largest (M > IO^^Mq) sys- 
tems at all redshifts, whereas we make slight under-predictions of 
the numbers of intermediate mass systems. Again, we attribute this 
success to our relaxation of the 'dissipationless collapse' scenario, 
whereby disc formation and mild star formation are assumed to oc- 
cur upon gaseous collapse at all epochs. However, we do also note 
that our under-prediction of the number density in the lower mass 
ranges is a cause for further analysis. 



3.3 Cosmological star formation rate density 

The cosmological star formation rate density (psfr) evolution of 
the Universe (i.e., the global rate of star formation as a function 
of redshift) is a key constraint for theoretical models of galaxy for- 
mation and cosmology, indicating a clear evolutionary link between 
the star forming properties of galaxy populations over different red- 
shifts. 

The 'Madau diagram', (Madau et al. 1996) has been used as 
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Figure 8. The evolution of the black hole mass function as compared to the determinations of Shankar et al. 2009. At z = we find encouraging agreements 
(left panel) across the observational range, at z = 1 we find that we slightly under-predict the number of high mass SMBHs (center panel) but note that this 
discrepancy may be related to the critical assumption of the validity of the Magorrian (1998) relation at 2 > 0. Finally we predict the mass function at 2 = 2 
(right panel), finding a further drop in normalisation and a lower cutoff. Yellow shaded regions represent the errors in the functional fits to observations (orange 
lines) and model eiTor bars represent the statistical uncertainties in mean averages, see Table Al. for tabulated data. 



a tool for constraining galaxy evolution models, however it's de- 
termination observationally is far from straightforward due to large 
systematic errors in extracting the SFR from luminosities, and cor- 
recting for dust obscuration and incompleteness. These factors led 
early determinations of the diagram to show a rapid increase by 
approximately an order of magnitude in psfr from < z < 1.5 
followed by a peak at z « 1.5 then a steady decline at higher red- 
shifts (see Madau et al. 1996 Fig.9). However, more sophisticated 
dust modeling and more complete samples al z > 1.5 resulted in 
revised estimates of the high-redshift decline; showing a relatively 
flat Psfr out to high-redshifts (Steidel et al. 1999). Combined with 
observations in the far-IR and sub-mm at intermediate and high 
redshifts (Hughes et al. 1998, Flores et al. 1999) the redshift depen- 
dence of the SFR density has become relatively well constrained to 

Originally theoretical attempts to predict the 'Madau diagram' 
were unable to model the correct evolution (see Cole et al. 1994), 
however, later works were able to match the results to a good ac- 
curacy (including the high-redshift decline, see Cole et al. 2000), 
and even within the latest generation of SAMs, a modest decline is 
observed between 2 > z > 5 (Somerville et al. 2008b), unlike the 
most recent observational constraints showing a near-flat evolution 
to 2 Ri 6 (see Hopkins, 2004). 

Within this work we utilize the observational compilation in 
Somerville et al. 2001, which discusses all the aforementioned sys- 
tematics and corrects for them accordingly (see references therein), 
also this work accounts for the standard cosmology. Shown in 
Fig.4, by fitting a cubic polynomial through the data shows (with 
considerable scatter), a general behavior of a rise in psfr from 
< 2 < 2 followed by a flattening at 2 < 2; < 3 and a 
slow decrease to higher redshifts, dropping to the 2 = value 
at 2 ~ 5.5. Outputting the total model Psfr, we see an overall 
agreement within the observational range, matching all the obser- 
vational features. We physically interpret the increase in psfr be- 
tween < 2 < 2 to several factors, increasingly rapid growth of 
DM haloes at higher redshifts allows more infalling material and at 
2 ~ 2 over the mass range of galactic haloes we have the syn- 



chronous formation of spheroid and disc components (since ap- 
proximately half of the haloes within our sample are in the 'fast 
collapse' phase and vice versa). Overall, above 2 ~ 2 we are dom- 
inated by the growth of spheroid-SMBH systems through the dis- 
sipative condensation of gas within DM haloes, therefore we typ- 
ically have higher psfr than predicted by SAMs constructed upon 
the disc-merger scenario. This naturally gives rise to a slow de- 
crease in Psfr to high redshifts. 

Also, we plot Psfr separated into both the bulge and disc com- 
ponents, since it has been suggested that they typical 'Hubble-type' 
morphological classification may be better understood as resulting 
from differing superpositions of spheroid and disc components (see 
Driver et al. 2006) where spheroids and discs form two separate 
classes each with their own distinct formation epochs and mecha- 
nisms. We find, in broad agreement with archeological studies, that 
Psfr is dominated by the spheroid component until 2 ~ 1, and be- 
comes progressively dominated by the disc component at 2 < 1, 
in accordance with the view that spheroids (and galaxy bulges) are 
typically 'red and dead', with old stellar populations, whereas discs 
show ongoing star formation over longer durations. 

3.4 Evolution of the galaxy scaling relations 

As a step beyond simply predicting the accumulation of stellar mat- 
ter within galaxies, theoretical models are able to make predictions 
about the dynamics and structure of galaxies which form, allow- 
ing for a further level of predictions and constraints from obser- 
vations. Initial observational advances in this direction came from 
studies of observable properties of galaxies. TuUy & Fisher, 1977, 
showed that a tight correlation existed between galaxy luminosity 
and maximum rotational velocity, these determinations have been 
confirmed and constrained in a number of latter works (see Haynes 
et al. 1999). 

The TuUy-Fisher relation (TFR) thus provides a link between 
luminous matter (stellar mass) and dynamical matter (total gravi- 
tational mass) of galaxies, providing strong constraints on the link 
between the underlying DM potential and the baryonic matter. Un- 
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fortunately however, theoretical attempts to interpret this relation 
within the framework of full SAMs have found many difficulties; 
offsets to within 30% are generally predicted by models (Cole et 
al. 2000) reinforcing the fact that simultaneous predictions of the 
stellar mass budgets and the TFR provide tight constraints on mod- 
els. This is further complicated since 'typical' SAMs make several 
assumptions and approximations in order to predict the maximum 
rotational velocity (see Cole et al. 2000). Within our model we di- 
rectly compute the rotation curve for the composite system given 
the density distributions of the DM halo, disc, bulge and central 
reservoir-SMBH system, providing us with a detailed output. Fol- 
lowing this, in accordance to observational methods, we output the 
value for the total rotation curve at 2.2 scale radii, typically corre- 
sponding to the 'peak' value for the rotation. Using this value, and 
plotting against the total galaxy mass in Fig. 5 we output the TFR's 
at three different redshifts. 

Comparing model results at z = to Bell & de Jong, 2001, 
and correcting for the stellar mass determinations due to different 
IMF choiceo we are able to make a good match within the obser- 
vational range (see also C09a), this indicates clearly that the model 
prescriptions which govern the baryon-to-DM ratios, and the struc- 
ture of DM haloes are producing the correct dynamical properties 
at fixed stellar mass. We note, however, that we have investigated 
the effects of DM halo contraction due to the condensation of bary- 
onic material, however, we find that we most accurately fit to ob- 
servational results without this effect (setting F = see Eqn.l3). 
This interesting finding has also been confirmed in several works 
focusing on the detailed structural properties of galaxies (see Dut- 
ton et al. 2007,2008), concluding that either a low M/L, or no halo 
contraction may be the only viable routes to achieving simultane- 
ous fits to the TFR and mass functions. Noting also the preliminary 
work presented in Mo & Mao, 2004, whereby haloes may become 
'pre-processed' due to an early rapid infall of matter, and ensu- 
ing mass outflow (through feedback) is able to considerably reduce 
halo concentrations (see §2.2 of Mo & Mao, 2004, see also El-Zant, 
Shlosman & Hoffman, 2001,2004, Tonini et al. 2006). We hope to 
further quantify these effects in a subsequent work. 

Outputting model results at higher redshifts we find little evo- 
lution between < 2; < 1, and, comparing results to the observa- 
tional determinations of Conselice et al. 2005 we find a good agree- 
ment, indicating that a general 'inside out' growth of discs, coupled 
with a growing DM halo results in galaxies that typically evolve 
along the TFR, thus showing little evolution. Despite the small ob- 
servational sample size (18 galaxies), we also show the z = 2 TFR 
as output by our model and compared to Cresci et al. 2009. Within 
their work, using the SINS survey (Forster Schreiber et al. 2006) 
which uses integral field spectroscopy to measure the dynamics of 
high-z galaxies, finding large rotating systems already in place at 
z > 2. We find that overall our model shows discrepancies with 
this data, however, we also note that during these epochs the stan- 
dard morphological sequence of galaxies is yet to be formed, and 
many galaxies within this sample are not in an equilibrium state but 
merely 'rotationally dominated' . 

Additionally to the TFR, we show in Fig.6 the evolution of the 
disc size relation. Basing on the work of Somerville et al. 2008a, 
who, motivated by detailed observations showing no significant 
evolution in the relationship between radial size and stellar mass 



^ Within Bell & de Jong, 2001, they take values for mass to light ratio 
which are approximately 30% lower than the Salpeter value, which we ac- 
count for when comparing stellar masses within this work. 
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Figure 9. The evolution of the gas-star fraction at different stellar masses, 
comparing model results to Baldry et al. 2008. By directly outputting the 
stellar and HI masses we are able to directly compare model results, find- 
ing good agreement at z = 0, we also show the evolution of this relation 
to z = 5 showing that low mass galaxies move onto the 2 = relation at 
early times, with the larger systems becoming more gas-iich at later times. 
The shaded yellow region represents the errors in the functional fit to obser- 
vations (orange lines) and error bars on model outputs representing Poisson 
uncertainties in the mean values, see Table A2. for tabulated data. 

from 2: « 1 to the present day, conducted a study of theoretical 
model predictions. By comparing our model predictions to results 
compiled from the Sloan Digital Sky Survey (SDSS), (Blanton et 
al. 2005, Somerville et al., 2008a) at « = we find a strong rela- 
tion between the disc scale radius (defined to be Rd = 0.5959i?Jfl 
see Courteau et al. 2007) and the disc stellar mass over the obser- 
vational range. At 2 = 1 however, we find that we systematically 
under predict the scale radius with an offset of ~ 1.5 compared 
to observations of GEMS galaxies (Rix et al., 2004, Caldwell et 
al., 2005, Somerville et al., 2008a), generating galaxies which are 
slightly too small at fixed stellar mass compared with observations. 
These subtle effects, however, are interesting since they appear to 
be showing that a scenario whereby an initial baryonic collapse and 
accompanying outflows thus lowering DM halo concentration may 
help to alleviate these issues (see Mo & Mao, 2004). 

3.5 Evolution of the black hole scaling relations 

Since the coevolution of a central SMBH and galaxy is a central 
importance mechanism in the formation and evolution of our sys- 
tem, we show here the evolution of the black hole scaling relations. 
In the local Universe it is now well established that most galactic 
nuclei host a central SMBH (Kormendy & Richstone, 1995). Each 
SMBH correlates strongly in mass with the global properties of the 
spheroid component of the host galaxy (Magorrian et al. 1998, Fer- 
rarese & Merritt, 2000), therefore a theoretical understanding of 
these relationships is of fundamental importance to galaxy forma- 
tion theories, and has been shown to provide the key in order to ac- 
count for the suppression of the formation of massive galaxies (see 

^ We note that there are sevearl cautionary details when converting half- 
light radii into stellar exponential scale lengths. Firstly, using sersic fitting 
methods may result in erroneous results for non-exponential discs, and sec- 
ondly, fitting using circular models generates size biases for inclined galax- 
ies. See Blanton et al., 2005 for more details. 
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Granato et al. 2004, Croton et al. 2006, Bower et al. 2006, Ciotti & 
Ostriker, 2007). In Fig.7 we compare model predictions of SMBH 
masses withi botli the velocity dispersion ((Tbujgefl and masses of 
the spheroid components. As can be seen, when comparing to the 
z = properties we find a close agreement to the observational 
constraints, both to the M-a relation (Tremaine et al. 2002) and 
the mass scaling (Haring & Rix, 2004). When comparing model 
predictions to higher redshifts, we utilize the constraints by Treu 
et al. 2007 & Woo et al. 2008 who used high resolution imaging 
in order to determine the SMBH masses and velocity dispersions 
at z — 0.36 and z = 0.57 respectively. Comparing these to model 
predictions we find a small offset, since we find that we do not have 
any significant evolution in the relations. However, in order to fur- 
ther constrain models, and determine whether the observed minor 
offset is physical or an artefact of increased scatter should allow us 
to further refine our computations (see Woo et al., 2008). Finally, 
for completeness, we show our prediction for the mass scalings at 
z — 1, again showing that there is no evolution in our models. 

Physically, this is due to the rapid growth of SMBH's within 
galaxies which fix the scaling relations, followed by periods of dor- 
mancy. We therefore expect scatter to increase at higher redshifts 
due to observations of galaxies which are still in the process of fix- 
ing their scaling relations. 

3.6 Black hole mass function evolution 

Utilising and exploiting the strong relationship between the central 
SMBH mass and the spheroid mass, along with methods in order to 
convert quasar (QSO) number counts into accreted mass densities 
onto central SMBH's, attempts to constrain a SMBH mass func- 
tion were initially conducted by Small & Blandford, 1992. Follow- 
ing this, several works ( see Salucci et al. 1999) related the lumi- 
nosity functions of local AGN's and of galaxy spheroids, in order 
to accurately determine the local mass function of SMBH's. Sev- 
eral further works using various techniques also made estimations 
(AUer & Richstone, 2002, Yu & Tremaine, 2002, McLure & Dun- 
lop, 2004, Marconi et al., 2004), however, several contrasting re- 
sults were produced. Shankar et al. 2004 later developed a robust 
method to determine the local SMBH mass function. In order to 
make estimations of the SMBH MF at 2 > we utilise the meth- 
ods outlined in Shankar, Bemardi & Haiman, 2009. In their work, 
they generate a SMBH mass function at different redshifts by map- 
ping the stellar mass function at z > onto a SMBH mass function 
through a Jacobian transformation, assuming the Magorrian (1998) 
relation remains valid at higher redshifts. 

In order to account for scatter, we convolve model outputs 
with a Gaussian scatter (0.3dex) and present the results in Fig. 8. 
As can be seen we match the z — mass function over ~ 4 orders 
of magnitude in mass, however at 2 = 1 we see model predic- 
tions fall slightly below the number density in the largest systems 
(PIbh > WH'Iq. Finally, we show our prediction for the SMBH 
mass function at 2 = 2. Assuming the Magorrian relation to be 
consistent with the local values is uncertain in these early epochs, 
and we advertise this as a direct prediction of the model. 

3.7 Stellar-to-Gas fraction evolution 

Outputting and analysing the evolution of cold gas within our 
model, we may make firm predictions to the total neutral gas frac- 

* We compute the velocity dispersion as Ubulge = O-GSVtiiiige 
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Figure 10. The evolution of the galaxy morphologies as a function of virial 
mass output by our model in the range < z < 5. Showing that at any 
epoch we have galaxies of all morphologies, but the characteristic transition 
mass decreases with increasing redshift. 

tion of galaxies as a function of stellar mass. As is shown in Fig.9, 
comparison to the compiled observations of Baldry, Glazebrook & 
Driver, 2008, show an agreement to observations over the entire 
mass range (4 magnitudes), this encouraging result highlights the 
accuracy in modeling the conversion of gaseous matter to stellar 
matter (see C09b for a detailed discussion). In order to make sev- 
eral predictions of our model, we show the evolution of this rela- 
tion at higher redshifts, finding that globally the gas fraction within 
galaxies increases at lower redshifts, with progressively lower mass 
galaxies being gas rich at higher redshifts. This is at variance with 
the general view, not yet observationally verified, that gas fractions 
increase with higher redshift. For instance, Stewart et al. 2009 as- 
sume such a positive z evolution for M gas /M stars, based the UV 
selected sample at z=2 by Erb et al 2006, and their estimate of the 
gas surface density, based on the Schimdt law. However this deter- 
mination is indirect and biased toward star forming gas rich sys- 
tems. We obtain different results from our model partly because in 
the prediction we include galaxies of all morphologies and evolu- 
tionary phase. In our model, an initial collapse causes a large in- 
flux of gas into haloes at high redshift, causing the formation of a 
spheroid-SMBH system within short timescales, with an effective- 
ness which increases with halo mass. Then high mass systems are 
soon stripped of gas due to QSO activity, leaving them dormant un- 
til a possible secondary disc growth at late times, whilst lower mass 
systems remain gas-rich because QSO feedback is incapable of re- 
moving gas. This results in high-z and high mass gas poor systems 
which are relatively dormant until conditions become sufficient to 
support the growth of a disc. As discs grow through the accretion 
of gas and relatively low star formation rates, the galaxies become 
progressively more gas-rich, resulting in the tight relation as ob- 
served locally. Future more direct and less biased determinations 
of gas fractions at 2 > will represent an interesting test for our 
model, whose predictions are tabulated in the Appendix. 



3.8 Morphological evolution 

Morphological classification of galaxies has been used as a pow- 
erful tool in order to separate galaxies into evolutionary cate- 
gories. Since Hubble (1926, 1936) defined the classic 'tuning fork' 
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diagram, little modification was required to achieve the modern 
scheme (see Sandage, 1961). It has become clear that morpholog- 
ical type defines more than merely the appearance of a galaxy, 
but also highlights general properties of the formation and evo- 
lution mechanisms which shape the final galaxy properties, it is 
understood that this is due to the fundamental underlying galaxy 
disc and spheroid components which superposed may generate the 
majority of morphological types (see Driver et al. 2006 for a de- 
tailed discussion). Conveniently, galaxy formation models typically 
separate galaxies into spheroid and disc components and then at- 
tempt to translate these into 'early' and 'late' type classifications 
through post-processed parameterisations, however, these are rela- 
tively subjective and may be somewhat arbitrarily chosen in order 
to 'filter' synthetic galaxy populations (see Cole et al. 2000). 

Still little is known about the evolution of galaxy morphology 
(see Parry, Eke & Frenk, 2009), however, it is clear that under a 
physically motivated galaxy formation scenario, whereby spheroids 
typically comprise of an old, single stellar population with little 
gas, and discs are a continuously evolving stellar population with a 
gas rich ISM, the bulge-disc ratio should indicate clearly with sev- 
eral morphological classifications of galaxies. With this in mind, 
we show in Fig. 10 the evolution of the bulge to total mass ratio 
output by our model, finding that, at 2 = we typically produce 
the observed relationship with low mass DM haloes hosting late- 
type, disc dominated galaxies, and high mass DM haloes hosting 
early-type, spheroid dominated galaxies, with a transition close to 
Li,, corresponding to ^ lO" - 10^^ Af©. We find that, at all red- 
shifts all morphologies are present, however, the transition between 
spheroid dominated galaxies drops to progressively lower masses. 
We hope to further investigate the morphological evolution, aiming 
to compare to bulge-disc decomposition studies at higher redshifts 
as future observations become available. 



3.9 Galactic Archeology 

Plotting the average stellar age of each galaxy against the z = Q 
mass of the galaxy and comparing to the results of Gallazzi et al. 
2005. 

Finally we compare the average stellar ages of galaxies of dif- 
ferent mass, the so called 'galactic archeology'. Observationally 
Gallazzi et al. 2005 used high resolution SDSS spectra in order 
to derive estimates for the ages and metallicities of ~ 170, 000 
galaxies through spectral and index fitting to a library of synthetic 
SEDs. These observations are shown in Fig. 1 1 , the large scatter be- 
ing attributed to the model-dependent age estimation of the stellar 
populations. As can clearly be seen, the mean stellar age of galax- 
ies decreases with decreasing mass, in contradiction with the naive 
'bottom-up' formation scenario, whereby we expect larger galax- 
ies to form at later times. This 'archeological downsizing' has been 
discussed in detail in Fontanot et al. 2009 (Fig.9), and we find, 
as with other current SAMs (Wang et al. 2008, Somerville et al. 
2008b, Monaco et al. 2007) that we are able to effectively predict 
the ages of the largest mass galaxies, but the low mass galaxies 
form and evolve too early, showing no clear signs of an archeo- 
logical downsizing. We therefore conclude that our model has dif- 
ficulties in predicting correctly the properties of the lowest mass 
galaxies, which typically form too early and are thus contain stellar 
populations which are too old. We view this as a significant lim- 
itation to our model (and to the aforementioned models) and this 
deserves further analysis. 
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Figure 11. The average age of the stellar masses of galaxies as a function 
of their 2 = stellar mass compared to the observations of Gallazzi et al. 
2005. We find, as with several other SAMs, that we are able to match the 
ai'cheological ages of the high mass galaxies, but low mass galaxies clearly 
form too early in our models, resulting in a significant over-prediction of the 
ages of these systems, e. The yellow shaded region represents the errors in 
the functional fit to observations (orange lines) and error bars in the model 
outputs represent the Poisson uncertainties in the mean averaged values. 



4 CONCLUSIONS 

Significant recent observational advances have allowed for un- 
precedented new constraints on the galaxy population, both locally 
and increasingly high redshifts. Motivated by the emerging phe- 
nomenological picture of galaxy evolution, we have presented a 
theoretical framework in order to interpret several observational 
constraints, finding a general agreement with several key results, 
some of which other SAMs find hard to reproduce. Within this 
work, we have expanded our two-phase galaxy formation model 
presented in (C09a, C09b) which constitutes a natural extension of 
the spheroid-SMBH coevolution model presented in Granato et al. 
2001, 2004 (see also Granato et al. 2006, Silva et al. 2005, Lapi 
et al. 2006, 2008, Mao et al. 2007). This model has been shown 
to naturally reproduce several key results, such as the properties 
of local elliptical galaxies, the sub-mm galaxy statistics, deep K- 
band survey results along with the local SMBH mass function and 
the statistics of high-z QSO's, the local gas fraction, HI and H2 
mass functions, stellar and baryonic mass functions, local lumi- 
nosity functions (separated into bulge and disc components), and 
the local Tully-Fisher relation. We inherit these results within this 
framework. 

The basic framework of the model presented here differs sig- 
nificantly from the typical 'disc-merger' scenario for galaxy forma- 
tion, challenging the assumption that gas cooling and condensation 
within DM haloes ubiquitously results in the dissipationless col- 
lapse onto a disc structure. By allowing for the direct infall onto a 
spheroid structure at early times we naturally generate large star- 
burst activity and thus are able to rapidly grow the largest galaxies 
at high-z, in agreement with many seemingly troublesome observa- 
tions. We note that the framework presented here is not incompat- 
ible with the 'standard' disc merger scenario, but we strongly rec- 
ommend the assumption of dissipationless collapse at any epoch 
to form disc structures should be further investigated within the 
latest SAMs (also noting that within hydrodynamic simulations an- 
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gular momentum dissipation is commonly seen, Zavala, Okamoto 
& Frank, 2008, Govemato et al., 2007). 

In order to include the most state-of-the-art processes thought 
relevant for galaxy formation, we include the effects of cold accre- 
tion flows, ionizing UV background radiation, a two-phase ISM, 
'radio mode' nuclear feedback, QSO and supernovae feedback, adi- 
abatic response of the DM halo to baryonic structure formation and 
disc stability criteria. This allows us to directly compare our re- 
sults with observations and other current SAMs. We find that, after 
accounting for the one-to-one relationship between DM halo and 
galaxy properties due to our simplified modeling, we are able to 
accurately reproduce the stellar mass function in the redshift range 
< z < 4 finding discrepancies only in the very lowest mass 
haloes (Fig. 2). We also show that we are able to match the evolution 
of massive galaxies from z — 5, attributing this to the early rapid 
growth allowed in our models through dissipative collapse onto a 
spheroid-SMBH structure (Fig.3), and we show the evolution of the 
cosmological SFR density (Fig.4), finding overall agreement with 
observations from 2; « 5, and also showing that SFR in spheroids 
dominates at z > 1, and at z < 1 the Universe favors quiescent 
disc growth. 

Focusing on the galaxy scaling relations, we show (in Fig.5) 
how the stellar mass Tully-Fisher relation shows little evolution 
from < 2; < 1, but then a marked difference at z = 2, and 
how the disc size evolution also shows Uttle evolution (Fig.6), this 
naturally results from our models since we naturally generate an 
'inside out' growth of discs, where the baryonic and DM evolve 
together. Also, since we place the mutual feedback between SFR 
and SMBH growth into central importance for the evolution of the 
most massive systems, we show in Fig.7 the evolution of the SMBH 
scaling relations, finding encouraging agreement at 2; = and little 
evolution to higher redshifts, also we show (Fig.8) the SMBH mass 
function compared to local estimates, and evaluate this at higher 
redshifts, tentatively comparing it to empirical fitting at z = 1 and 
predicting the evolution to « = 2. In order to further highlight the 
evolutionary differences in our model, we show in Fig.9 the evo- 
lution of the stellar-gas fraction, showing how the 2 = relation 
is constructed, and predict the growth of morphological types in 
Fig. 10. Finally, motivated by recent determinations of mean stel- 
lar ages of galaxies, 'archeological downsizing' has been noted in 
the literature, in Fig. 11 we show the mean stellar ages of galaxies 
as a function of their 2 = stellar mass, finding, as with other 
SAMs, that the smallest systems form too early in our framework 
(see Fontanot et al. 2009), we confirm that at present this is a robust 
challenge to all current models. 

In summary, under our proposed framework we are able to si- 
multaneously reproduce the vast majority of key observational con- 
straints on galaxy formation in the range < 2; < 5, concentrating 
particularly on several plots which are notoriously troublesome for 
SAMs to reproduce, finding minor discrepancies between model 
and observation mainly where observational results are not con- 
straining and subject to large potential biases. We therefore regard 
this as a large success of our model. Coupled with the successes of 
this framework in previous papers, focusing on the detailed proper- 
ties of galaxies both locally and at high-z, and further advances in 
observational constraints (particularly with resolved spectroscopy) 
we hope to further constrain the detailed processes governing the 
evolution of baryonic matter within evolving DM haloes. 

We also note that there are several points of tension within 
our model, manifesting within the lowest mass systems. At least 
in part, these tensions are Ukely due to oversimplified star forma- 
tion recipes. Indeed, as discussed in detail in C09b, by including 



a two-phase ISM, evoking a SFR related to the molecular gas sur- 
face density, and including ionizing UV background suppression 
we are able to significantly reduce the number of low mass satellite 
galaxies. 

However, the main limitation to our simplified computations 
is that we neglect the environmental effects (tidal stripping and 
harassment) due to our single MAH approach. We note that the 
lowest mass galaxies within our observational range are typically 
the ones most likely to be embedded within larger structures and 
thus will be subject to external effects (see Mo et al., 2005), by ac- 
covmting for these we hope to achieve closer matches to the faint- 
end slope of the stellar mass functions and this may also help to 
alleviate model discrepancies with the 'archeological downsizing' 
problem, by pre-heating material and preventing it infalling on the 
lowest mass haloes at early times due to embedding within larger 
structures. We thus hope to explore these effects within a subse- 
quent work, noting again that our prescriptions are not mutually 
incompatible with current SAMs, but simply require modification 
to loosen the assumption of ubiquitous dissipationless gaseous col- 
lapse and naturally ease several tensions between theory and obser- 
vation. 

We therefore advocate the exploration of the evolution of an- 
gular momentum in simulations, and hope to conduct a further 
analysis into how this may be physically modeled within a self- 
consistent semi-analytical model, basing on the successes of this 
simpUfied approach. 



APPENDIX A: TABULATED MODEL PREDICTIONS 

Within this appendix we tabulate the results for the black hole mass 
function at 2 = 2 and the stellar-to-gas fraction up to 2; = 5 where 
currently there are no strong constraints, but with future studies, 
these model predictions may be compared with observations. 
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